System and method for generating mr myocardial perfusion maps without user interaction

ABSTRACT

A method for automatically generating a myocardial perfusion map from a sequence of magnetic resonance (MR) images includes determining a region of interest (ROI) in a reference frame selected from a time series of myocardial perfusion MR image slices, registering each image slice in the time series of slices to the reference frame to obtain a series of registered ROIs, and using the series of registered ROIs to segment endo- and epi-cardial boundaries of a myocardium in the ROI.

CROSS REFERENCE TO RELATED UNITED STATES APPLICATIONS

This application claims priority from “Image Analysis Algorithms forGenerating MR Myocardial Perfusion Maps without User Interaction”, U.S.Provisional Application No. 61/055,180 of Sun, et al, filed May 22,2008, the contents of which are herein incorporated by reference intheir entirety.

TECHNICAL FIELD

This disclosure is directed to the automatic generation of magneticresonance (MR) myocardial perfusion maps.

DISCUSSION OF THE RELATED ART

In physiology, perfusion is the process of nutritive delivery ofarterial blood to a capillary bed in biological tissue. Perfusionscanning is the process by which this perfusion can be observed,recorded and quantified. Being able to observe and quantify perfusion inthe human body has been an invaluable step forward in medicine. With theability to ascertain data on the blood flow to vital organs such as theheart and the brain, physicians can make quicker and more accuratechoices on treatment for patients. Perfusion scanning encompasses a widerange of medical imaging modalities, such as computed tomography (CT)and magnetic resonance imaging (MRI). There are several differenttechniques of detecting perfusion parameters with the use of MRI, themost common being dynamic susceptibility contrast imaging (DSC-MRI) andarterial spin labeling (ASL).

Dynamic perfusion magnetic resonance imaging (MRI) has demonstratedgreat potential for diagnosing cardiovascular and renovascular diseases.In cardiac perfusion, the contrast agent passes through the rightventricle to the left ventricle and then perfuses into the myocardium.In dynamic perfusion MRI, the organ under study is scanned rapidly andrepeatedly following a bolus injection of a contrast agent. Thekinematics of the contrast agent are reflected in the intensity changesof the obtained time series of MR images. Analysis of the dynamicbehavior of the signal intensity can provide valuable functionalinformation.

Myocardial perfusion scanning is a procedure that evaluates many heartconditions from coronary artery disease to hypertrophic cardiomyopathy.The function of the myocardium is also evaluated by calculating the leftventricular ejection fraction of the heart. Myocardial perfusion MRimages are generally compromised by low signal-to-noise ratio (SNR) andpoor spatial resolution relative to the size of potential ischemicregions due to rapid imaging on the order of 100-150 ms per image. Forquantitative or semi-quantitative analysis of these images, parametricmaps are often generated of parameters, such as upslope, peak intensity,and area under the curve. Current methods for parametric map generationrequire significant user interaction for registration or segmentationsteps.

SUMMARY OF THE INVENTION

Exemplary embodiments of the invention as described herein generallyinclude methods and systems for automatically generating myocardialperfusion maps from MRI volumetric datasets, by determining of a regionof interest (ROI) around the heart, registering the ROI across an MRvolumetric image sequence, and segmenting the ROI to detect the endo-and epicardial boundaries of the myocardium.

According to an aspect of the invention, there is provided a method forautomatically generating a myocardial perfusion map from a sequence ofmagnetic resonance (MR) images, including determining a region ofinterest (ROI) in a reference frame selected from a time series ofmyocardial perfusion MR image slices, registering each image slice inthe time series of slices to the reference frame to obtain a series ofregistered ROIs, and using the series of registered ROIs to segmentendo- and epi-cardial boundaries of a myocardium in the ROI.

According to a further aspect of the invention, determining a region ofinterest (ROI) in a reference frame comprises cropping slices in theseries of slices according to a temporal intensity variation at eachpixel to obtain a global mask containing pixels exhibiting significanttemporal intensity variation, selecting a range of candidate slices withsuperior contrast, detecting the LV blood pool and forming an LV mask inall candidate slices, selecting one candidate slice with a detected LVblood pool as the reference frame, detecting the right ventricle (RV) inthe reference frame and forming an RV mask, estimating a thickness ofthe myocardium based on the detected LV and RV masks, and defining theROI as a bounding box containing the LV, RV, and myocardium.

According to a further aspect of the invention, detecting the LV bloodpool comprises, for each candidate slice, selecting an intensitythreshold for the blood, using the threshold to obtain a binary imageand removing small connected components from the binary image.

According to a further aspect of the invention, selecting one candidateslice as the reference frame comprises computing shape features for eachremaining connected component to select a least eccentric, mostcircular, and most convex connected component, where the selectedconnected component among all among all candidate slices corresponds tothe LV blood pool.

According to a further aspect of the invention, the blood intensitythreshold is selected through Gaussian clustering.

According to a further aspect of the invention, the shape featuresinclude an eccentricity computed as a ratio between minimum and maximumradii on the connected components, a convexity computed as a ratiobetween an area of the connected component and an area of its convexhull, and a circularity of the convex hull defined as 4π times the areaof the convex hull over a square of a perimeter of the convex hull.

According to a further aspect of the invention, when two connectedcomponents have similar shape features, a component with a larger areais selected.

According to a further aspect of the invention, detecting the RVcomprises finding a connected component in the binary image that isclose to the LV and has a large area.

According to a further aspect of the invention, segmenting endo- andepi-cardial boundaries of a myocardium in the ROI comprises initializingthe endo-cardium and epi-cardium boundaries from the LV mask andmyocardial thickness estimate, detecting the endo-cardium by minimizingan energy functional, masking out a left-ventricle (LV) blood poolinside of the endo-cardium from the image series, and detecting theepi-cardium by minimizing another energy functional, where each energyfunctional incorporates correlation information among pixels in a sameslice as well as temporal correlations across the slices in the timeseries, prior shape information and anatomical constraints.

According to a further aspect of the invention, the energy functionalfor the endo-cardium isE(C₁)=E_(seq)(C₁)+λE_(img)(C₁)+μE_(int)(C₁)+κE_(shape)(C₁), where C₁ isthe endo-cardial boundary, λ, μ, and κ are parameters that controlweights of different terms, where

${{E_{seq}(C)} = {{\int_{\Omega_{i\; n}}{{w\left( {x,y} \right)}{\begin{bmatrix}\begin{matrix}{{\alpha \; {{dis}^{2}\left( {{I\left( {x,y} \right)},{\overset{\_}{I}}_{i\; n}} \right)}} +} \\{{\beta \; {{dis}^{2}\left( {{I\left( {x,y} \right)},{\overset{\_}{I}}_{prior}} \right)}} +}\end{matrix} \\{\gamma \; {{dis}^{2}\left( {{\overset{\_}{I}}_{i\; n},{\overset{\_}{I}}_{prior}} \right)}}\end{bmatrix}}{x}{y}}} + {\int_{\Omega_{out}}{\left\lbrack {{w\left( {x,y} \right)}{{dis}^{2}\left( {{I\left( {x,y} \right)},{\overset{\_}{I}}_{out}} \right)}} \right\rbrack {x}{y}}}}},$

where

${{w\left( {x,y} \right)} = \frac{\sum\limits_{k = 1}^{K}\left\lbrack {{I_{k}\left( {x,y} \right)} - {I_{base}\left( {x,y} \right)}} \right\rbrack}{\max_{({x,y})}{{\sum\limits_{k = 1}^{K}\left\lbrack {{I_{k}\left( {x,y} \right)} - {I_{base}\left( {x,y} \right)}} \right\rbrack}}}},$

I_(k)(x,y) is the intensity of a pixel (x,y) in a kth slice of the timeseries, I_(max)(x, y) is a baseline intensity for the pixel (x, y), I(x,y)=[I₁(x, y), I₂(x, y), . . . , I_(K)(x, y)]^(T) is a intensity-timecurve for each pixel (x, y), Ī_(in) and Ī_(out) are averageintensity-time vectors for region Ω_(in) and Ω_(out), respectively,Ω_(in) and Ω_(out) are the regions inside and outside boundary C₁,respectively, Ī_(prior) is a prior known average intensity-time curveinside boundary C₁,

$\begin{matrix}{{{dis}^{2}\left( {I_{1},I_{2}} \right)} = {\frac{1}{2} - \frac{\langle{{\overset{\sim}{I}}_{1},{\overset{\sim}{I}}_{2}}\rangle}{2{{\overset{\sim}{I}}_{1}}_{2}{{\overset{\sim}{I}}_{2}}_{2}}}} \\{= \frac{{{{\overset{\sim}{I}}_{1}}_{2}{{\overset{\sim}{I}}_{2}}_{2}} - {\langle{{\overset{\sim}{I}}_{1},{\overset{\sim}{I}}_{2}}\rangle}}{2{{\overset{\sim}{I}}_{1}}_{2}{{\overset{\sim}{I}}_{2}}_{2}}}\end{matrix}$

where Ĩ is a mean removed vector, and α, β, and γ are positive scalarsthat lie in [0, 1] and satisfy α+β+γ=1,

${{E_{img}(C)} = {{\int_{\Omega_{i\; n}}{\begin{bmatrix}\begin{matrix}{{\alpha \left( {{I_{r}\left( {x,y} \right)} - m_{i\; n}} \right)}^{2} +} \\{{\beta \left( {{I_{r}\left( {x,y} \right)} - m_{prior}} \right)}^{2} +}\end{matrix} \\{\gamma \left( {m_{i\; n} - m_{prior}} \right)}^{2}\end{bmatrix}{x}{y}}} + {\int_{\Omega_{out}}{\begin{pmatrix}{{I_{r}\left( {x,y} \right)} -} \\m_{out}\end{pmatrix}^{2}{x}{y}}}}},$

where I_(r) is the reference frame, m_(in) and m_(out) are averageintensities inside Ω_(in) and Ω_(out), respectively and m_(prior) is apriori known average intensity inside curve C₁, E_(int)=∫_(C) ₁ ds, anarc length of boundary C₁, and E_(shape)=∫_(C) ₁ D²(x, y)ds, where D(x,y) is a distance between a pixel (x, y) in an elliptical prior shape andthe estimated boundary C₁.

According to a further aspect of the invention, the method includesweighting each pixel on the endo-cardial boundary according to itsdistance to a convex hull of endo-cardial boundary, where pixels closerto the convex hull have a greater weight.

According to a further aspect of the invention, the energy functionalfor the epi-cardium isE(C₂)=E_(seq)(C₂)+λE_(img)(C₂)+μE_(int)(C₂)+ηE_(shape)(C₂)+E_(con)(C₂),where C₂ is the epi-cardial boundary, λ, μ, and κ are parameters thatcontrol weights of different terms, where

${{E_{seq}(C)} = {{\int_{\Omega_{{i\; n}\;}}{{{w\left( {x,y} \right)}\begin{bmatrix}\begin{matrix}{{\alpha \; {{dis}^{2}\left( {{I\left( {x,y} \right)},{\overset{\_}{I}}_{i\; n}} \right)}} +} \\{{\beta \; {{dis}^{2}\left( {{I\left( {x,y} \right)},{\overset{\_}{I}}_{prior}} \right)}} +}\end{matrix} \\{\gamma \; {{dis}^{2}\left( {{\overset{\_}{I}}_{i\; n},{\overset{\_}{I}}_{prior}} \right)}}\end{bmatrix}}{x}{y}}} + {\int_{\Omega_{out}}{\left\lbrack {{w\left( {x,y} \right)}{{dis}^{2}\left( {{I\left( {x,y} \right)},{\overset{\_}{I}}_{out}} \right)}} \right\rbrack {x}{y}}}}},$

where

${{w\left( {x,y} \right)} = \frac{\sum\limits_{k = 1}^{K}\left\lbrack {{I_{k}\left( {x,y} \right)} - {I_{base}\left( {x,y} \right)}} \right\rbrack}{m\; {ax}_{({x,y})}{{\sum\limits_{k = 1}^{K}\left\lbrack {{I_{k}\left( {x,y} \right)} - {I_{base}\left( {x,y} \right)}} \right\rbrack}}}},$

I_(k)(x, y) is the intensity of a pixel (x, y) in a kth slice of thetime series, I_(base)(x, y) is a baseline intensity for the pixel (x,y), I(x, y)=[I₁(x, y), I₂(x, y), . . . , I_(K)(x, y)]^(T) is aintensity-time curve for each pixel (x, y), Ī_(in) and Ī_(out) areaverage intensity-time vectors for region Ω_(in) and Ω_(out),respectively, Ω_(in) and Ω_(out) are regions inside and outside boundaryC₂, respectively, and Ω_(in) excludes a region inside the endo-cardiumboundary C₁, Ī_(prior) is a prior known average intensity-time curveinside boundary C₂,

$\begin{matrix}{{{dis}^{2}\left( {I_{1},I_{2}} \right)} = {\frac{1}{2} - \frac{\langle{{\overset{\sim}{I}}_{1},{\overset{\sim}{I}}_{2}}\rangle}{2{{\overset{\sim}{I}}_{1}}_{2}{{\overset{\sim}{I}}_{2}}_{2}}}} \\{= \frac{{{{\overset{\sim}{I}}_{1}}_{2}{{\overset{\sim}{I}}_{2}}_{2}} - {\langle{{\overset{\sim}{I}}_{1},{\overset{\sim}{I}}_{2}}\rangle}}{2{{\overset{\sim}{I}}_{1}}_{2}{{\overset{\sim}{I}}_{2}}_{2}}}\end{matrix}$

where Ĩ is a mean removed vector, and α, β, and γ are positive scalarsthat lie in [0, 1] and satisfy α+β+γ=1,

${{E_{img}(C)} = {{\int_{\Omega_{i\; n}}{\begin{bmatrix}\begin{matrix}{{\alpha \left( {{I_{r}\left( {x,y} \right)} - m_{i\; n}} \right)}^{2} +} \\{{\beta \left( {{I_{r}\left( {x,y} \right)} - m_{prior}} \right)}^{2} +}\end{matrix} \\{\gamma \left( {m_{i\; n} - m_{prior}} \right)}^{2}\end{bmatrix}{x}{y}}} + {\int_{\Omega_{out}}{\left( {{I_{r}\left( {x,y} \right)} - m_{out}} \right)^{2}{x}{y}}}}},$

where I_(r) is the reference frame, m_(in) and m_(out) are averageintensities inside Ω_(in) and Ω_(out), respectively and m_(prior) is apriori known average intensity inside curve C₁, E_(int)=∫_(C) ₂ ds, anarc length of boundary C₂, E_(shape)=∫_(C) ₂ D²(x, y)ds, where D(x, y)is a distance between pixel (x, y) and the estimated boundary C₂, andE_(con)=∫_(C) ₂ f[d₁(x, y)]ds, where d₁(x, y) is a distance between apoint (x, y) on C₂ and C_(1,) and

${f\left\lbrack {d_{1}\left( {x,y} \right)} \right\rbrack} = \left\{ {\begin{matrix}{0,{{{if}\mspace{14mu} d_{m\; i\; n}} \leq {d_{1}\left( {x,y} \right)} \leq d_{{ma}\; x}}} \\{{\infty,{otherwise}}\mspace{149mu}}\end{matrix},} \right.$

where d_(min) and d_(max) are minimum and maximum allowed myocardiumthicknesses, respectively.

According to another aspect of the invention, there is provided a methodof automatically generating a myocardial perfusion map from a sequenceof magnetic resonance (MR) volumetric images, including providing a timeseries of registered myocardial perfusion MR image slices havingreference image with a masked left ventricle (LV), detecting the rightventricle (RV) in the reference frame and estimating a thickness of themyocardium based on the detected LV and RV, initializing theendo-cardium and epi-cardium boundaries from the LV mask and myocardialthickness estimate, segmenting the endo-cardium by minimizing an energyfunctional, masking out a left-ventricle (LV) blood pool inside of theendo-cardium from the image series, and segmenting the epi-cardium byminimizing another energy functional, where each energy functionalincorporates correlation information among pixels in a same slice aswell as temporal correlations across the slices in the time series,prior shape information and anatomical constraints.

According to a further aspect of the invention, providing a time seriesof registered myocardial perfusion MR image slices comprises providing atime series of myocardial perfusion MR image slices, determining aregion of interest (ROI) in a reference frame selected from the timeseries, and registering each image slice in the time series of slices tothe reference frame to obtain a series of registered ROIs, where theseries of registered ROIs are used to segment endo- and epi-cardialboundaries of a myocardium in the ROI.

According to another aspect of the invention, there is provided aprogram storage device readable by a computer, tangibly embodying aprogram of instructions executable by the computer to perform the methodsteps for automatically generating a myocardial perfusion map from asequence of magnetic resonance (MR) images.

BRIEF DESCRIPTION OF THE DRAWINGS

FIG. 1 is a flowchart of a method for mapping myocardial perfusion MRimages, according to an embodiment of the invention.

FIG. 2 is a flowchart of a method for determining a region-of-interestin a reference frame, according to an embodiment of the invention.

FIG. 3 depicts papillary muscles in the heart, according to anembodiment of the invention.

FIGS. 4( a)-(b) illustrate the intensity-time curves of 4 pixels acrossan image sequence, according to an embodiment of the invention.

FIG. 5 illustrates the blood pool, myocardium, and background regions,according to an embodiment of the invention.

FIG. 6 illustrates sample results of a region-of-interest (ROI)determination, according to an embodiment of the invention.

FIG. 7 illustrates sample segmentation results, according to anembodiment of the invention.

FIG. 8 is a block diagram of an exemplary computer system forimplementing a method for automatically generating magnetic resonance(MR) myocardial perfusion maps, according to an embodiment of theinvention.

DETAILED DESCRIPTION OF EXEMPLARY EMBODIMENTS

Exemplary embodiments of the invention as described herein generallyinclude systems and methods for automatically generating magneticresonance (MR) myocardial perfusion maps. Accordingly, while theinvention is susceptible to various modifications and alternative forms,specific embodiments thereof are shown by way of example in the drawingsand will herein be described in detail. It should be understood,however, that there is no intent to limit the invention to theparticular forms disclosed, but on the contrary, the invention is tocover all modifications, equivalents, and alternatives falling withinthe spirit and scope of the invention.

As used herein, the term “image” refers to multi-dimensional datacomposed of discrete image elements (e.g., pixels for 2-D images andvoxels for 3-D images). The image may be, for example, a medical imageof a subject collected by computer tomography, magnetic resonanceimaging, ultrasound, or any other medical imaging system known to one ofskill in the art. The image may also be provided from non-medicalcontexts, such as, for example, remote sensing systems, electronmicroscopy, etc. Although an image can be thought of as a function fromR³ to R, the methods of the inventions are not limited to such images,and can be applied to images of any dimension, e.g., a 2-D picture or a3-D volume. For a 2- or 3-dimensional image, the domain of the image istypically a 2- or 3-dimensional rectangular array, wherein each pixel orvoxel can be addressed with reference to a set of 2 or 3 mutuallyorthogonal axes. The terms “digital” and “digitized” as used herein willrefer to images or volumes, as appropriate, in a digital or digitizedformat acquired via a digital acquisition system or via conversion froman analog image.

According to an embodiment of the invention, the following steps in theanalysis of myocardial perfusion MR images can be automated. These stepsare depicted in the flowchart of FIG. 1. Referring to the figure, givena time series of myocardial perfusion MR images, an automated methodbegins at step 11 by determining a region of interest (ROI), i.e., abounding box around the heart. At step 12, registration of the ROIacross an image sequence is performed, and, at step 13, the ROI issegmented to detect the endo- and epicardial boundaries of themyocardium.

Determination of ROI

Given a time-series of myocardial perfusion MR images of a single slice,it is desired to automatically select a reference frame and generate arectangular ROI that tightly contains the left ventricle in thereference frame. The ROI is then used for registration across framesfollowed by the segmentation of the endo- and epicardial contours. Anexemplary, non-limiting time series would comprise about 40-60 MRframes.

This task is challenging for several reasons. The contrast of the MRimages varies significantly, as images acquired before the contrastagent washes in have poor contrast, while images acquired after thewash-in have much better contrast. Hence a selected reference frameshould have good contrast. In addition, the left ventricle blood pooland myocardium are not the only structures whose intensities change overtime, as the right ventricle and other structures also respond to thecontrast agent. Finally, the image sequence often experiencessignificant patient breathing motion. As a result, techniques such asmaximum intensity projection to produce a feature image for thedetection of the left ventricle (LV) cannot be used.

Among all the images in the time-series, there generally exists a fewimages that have good contrast. In these images, the intensity of theleft ventricle blood pool is higher than the rest of the structures. Inaddition, the shape of the left ventricle blood pool is alwaysellipse-like if not circular. For these reasons, according to anembodiment of the invention, the heart can be located by finding theleft ventricle (LV) blood pool in a high contrast images selected as areference image.

FIG. 2 is a flowchart of a method for determining a region-of-interestin a reference frame, according to an embodiment of the invention.Referring now to the figure, given an MR image sequence, a method beginsat step 21 by cropping images according to the temporal intensityvariation at each pixel to obtain a global mask containing pixelsexhibiting significant temporal intensity variation. According to theaverage intensity within the global mask at each frame, a range offrames with good contrast is selected at step 22. At step 23, the LVblood pool is detected in all candidate frames. Detecting the blood poolbegins by, for each selected frame, selecting an intensity threshold forthe blood through Gaussian clustering, and obtaining a binary image bythresholding. After removing small connected components, three shapefeatures are computed for each remaining connected component: (1) theeccentricity is computed as the ratio between the minimum and themaximum radii on the connected component; (2) the convexity is computedas the ratio between the area of the connected component and the area ofits convex hull; and (3) the circularity of the convex hull is definedas 4π times its area over the square of its perimeter. The leasteccentric, most circular, and most convex connected component is sought.The winning connected component among all frames is assumed tocorrespond to the LV blood pool. In addition, when two connectedcomponents have similar shape features, the component with a larger areais favored to avoid mistaking the aorta as the heart. At step 24, theframe containing the winning component is selected as the referenceframe, and the reference frame is used to estimate the thickness of themyocardium. The right ventricle (RV) is then found at step 25. Once theLV has been detected, the RV is then the connected component in thebinary image that is close to the LV and at the same time has a fairlylarge area. The thickness of the myocardium can be estimated based onthe detected LV and RV masks. An appropriate margin is set according tothe estimated thickness as well as the area of the LV. Once the masksfor the LV and RV blood pool and the estimated thickness of themyocardium have been determined, one can determine at step 26 a boundingbox that defines an ROI containing the LV.

Registration of the ROI

For quantitative or semi-quantitative analysis of cardiac perfusionimages, registration should be performed to minimize the errors due torespiratory motion. According to an embodiment of the invention, a rigidbody registration of the ROI is performed for each slice position. Eachindividual time frame is registered with respect to a reference image, aprocess sometimes referred to as tracking the ROI.

An exemplary registration method according to an embodiment of theinvention is as follows. This algorithm is disclosed in the assignee'spending application, “Contrast-invariant registration of cardiac andrenal magnetic resonance perfusion images”, U.S. patent application Ser.No. 11/078,035, filed Mar. 11, 2005, the contents of which are hereinincorporated by reference in their entirety. An algorithm begins, givenan ROI in one frame, by finding the best match to the selected ROI inthe other frames. According to an embodiment of the invention, it isassumed that the motion is reduced to translation with integer pixelshifts. Registering the ROI is achieved by simple template matching.Since the orientations of the edges along tissue boundaries are alwaysparallel across the image sequence, despite the fact that the relativeintensities between tissues vary with time, a template defined by theorientation of the image gradient is chosen.

In this formulation, the image on which the ROI is manually cropped iscalled the reference frame. Let θ_(r)(x, y) and M_(r)(x, y) representthe direction and the magnitude, respectively, of the image gradient atpixel (x, y) in the reference frame. One can obtain θ_(r)(x, y) andM_(r)(x, y) using, for example, a Sobel edge detector as is known in theart, although any edge detector may be used according to otherembodiments of the invention. The set of pixels in the ROI is denoted asR={(x, y)|x_(a)≦x≦x_(b), y_(a)≦y≦y_(b)}, where (x_(a), y_(a)) and(x_(b), x_(b)) are the diagonal points that specify the bounding box ofthe ROI. Let θ₁(x, y) denote the edge orientation and M_(c)(x, y) theedge magnitude at pixel (x, y) in the current frame. For each offsetpair (dx, dy), the angle difference Δθ(x, y; dx, dy) and a weightfunction w(x, y; dr, dy) are defined by:

$\begin{matrix}{{{\Delta \; {\theta \left( {x,{y;{x}},{y}} \right)}} = {{\theta_{c}\left( {{x + {x}},{y + {y}}} \right)} - {\theta_{r}\left( {x,y} \right)}}},} & (1) \\{{w\left( {x,{y;{x}},{y}} \right)} = {\frac{{M_{c}\left( {{x + {x}},{y + {y}}} \right)}{M_{r}\left( {x,y} \right)}}{\sum\limits_{{({x,y})} \in R}{{M_{c}\left( {{x + {x}},{y + {y}}} \right)}{M_{r}\left( {x,y} \right)}}}.}} & (2)\end{matrix}$

A similarity metric is introduced:

$\begin{matrix}{{S\left( {{x},{y}} \right)} = {\sum\limits_{{({x,y})} \in R}{{w\left( {x,{y;{x}},{y}} \right)}{{\cos \left( {2\Delta \; {\theta \left( {x,{y;{x}},{y}} \right)}} \right)}.}}}} & (3)\end{matrix}$

Note that S(dx, dy) is the weighted average of the values of cos(2Δθ)over the ROI, and its value lies in the interval of [−1, 1]. The cosineof the double angle is used to handle contrast inversions that commonlyoccur in perfusion studies. In addition, the weight function is chosenas the normalized product of the edge magnitudes because it is desirablefor the ROI to be attracted to strong edges in the current frame thatare mostly parallel to the strong edges in the reference frame. Theexemplary similarity metric is invariant to rapidly changing imagecontrast, in the sense that its value is insensitive to changes in thecontrast as long as the edge orientations are nearly parallel to thoseof the template. The integer shifts (dx*, dy*) that maximize S aredetermined by exploring all possible solutions (dx, dy) over areasonable search space. Note that the value of S(dx*, dy*) also plays arole as a confidence measure. To improve the robustness of thealgorithm, both the previous frame and the reference frame are used astemplates, choosing the match with maximum similarity metric. To furtherimprove the accuracy of the registration results, a local affinetransformation of the heart is estimated by incorporating the knowledgeof the contours delineating organ boundaries in the reference frame.

Segmentation of the ROI

The registered ROI sequence may then used to detect the inner and outerboundaries of the myocardium. However, in dynamic MR perfusion imagesequences, it can be challenging to distinguish different anatomicalstructures based on each individual frame, mostly due to a lack ofcontrast along some boundary segments. However, successful registrationof the ROI enables the use of the entire MR image sequence in thesegmentation. To use all the information available, an imagesegmentation method according to an embodiment of the invention uses notonly the spatial information provided by each single image in the MRimage sequence but also the temporal information available, inparticular, the intensity-time curves. This is motivated by thefollowing observation.

FIGS. 4( a)-(b) illustrates the intensity-time curves of 4 pixels 41,42, 43, 44 across the image sequence respectively located at the rightventricle 41, left ventricle 43, myocardium 42, and background 44. FIG.4( a) depicts the selected pixels, while FIG. 4(b) graphs theintensity-time curves for the selected pixels. These dynamic sequencesare all quite distinct. For example, the intensity of the rightventricle and the left ventricle pixels show pronounced oscillations andexhibit different wash-in timings, whereas the intensity of thebackground pixel is very noise like. This suggests that one candistinguish between the heart and other anatomical structures based ontemporally distinct patterns of their intensity-time curves. Inaddition, structures within the heart, such as the left ventricle cavityand myocardium can also be distinguished based on their temporalresponses.

A segmentation algorithm according to an embodiment of the invention canexploit all the information available, spatially or intra-image, as wellas temporally or across the images of the sequence, and can be tuned tothe specificity of cardiac MR perfusion images. Since the contrast atthe endocardial boundary is generally better than that at the epicardialboundary, an algorithm according to an embodiment of the invention firstsegments the endocardium and then segments the epicardium utilizing theresults of endocardium segmentation.

A segmentation algorithm according to an embodiment of the inventionincludes 2 steps: first, detecting the inner boundary (endo); then,masking out the LV blood pool (inside of the inner boundary) from theimage sequence and detecting the epi-cardium. Contours are initializedusing the LV blood pool and myocardial thickness estimates obtainedabove. According to an embodiment of the invention, segmentation isformulated as an energy minimization problem. Since the intensity-timeprofiles for pixels in the LV, RV and myocardium are quite different,energy functionals are introduced that exploit the difference in thedynamics of the temporal signals associated with distinct pixels. Theenergy functionals for the two segmentation tasks follow a general form,but with special variations in either case to address their uniquechallenges. An energy-based image segmentation algorithm according to anembodiment of the invention uses the correlation information amongpixels in the same image as well as the temporal correlation across theimages in the sequence. It also exploits prior shape information andanatomical constraints.

Because the two boundary contours are initialized according to the maskfor the LV blood pool and the estimated thickness of the myocardium, theinitial contours are very close to the true boundaries. Therefore, onlya small number of iterations are needed for the energy minimization stepto converge.

NOTATION: Let I_(k): ΩεR²→R denote the k^(th) frame in the imagesequence of the ROI to be segmented, k=1, 2, . . . , K, and I_(k)(x, y)represent the intensity of the pixel at (x, y). Here K represents thetotal number of the frames in the image sequence. Let C ε Q be a closedcontour. Let C₁ and C₂ denote respectively the evolving endocardial andepicardial contours in the reference frame. These two contours partitionthe image domain Q into three disjoint regions that correspondrespectively to the blood pool, the myocardium, and the background. FIG.5 illustrates the blood pool Ω_(a)=inside(C₁), the myocardiumΩ_(b)=outside(C₁) ∩ inside(C₂), and the background regionsΩ_(c)=outside(C₂).

ENERGY FUNCTIONAL: A general energy functional according to anembodiment of the invention is defined as:

E(C)=E _(seq)(C)+λE _(img)(C)+μE _(int)(C)+κE _(shape)(C),   (4)

where E_(seq) accounts for the temporal dynamics in the images sequence,E_(img) incorporates the spatial information in the reference frame,E_(int) is the internal energy that regulates the smoothness of thecontours, E_(shape) penalizes the shape difference between the evolvingcontour and a prior shape model, and λ, μ, and κ are parameters thatcontrol the weights of different terms. Next, each of these energy termswill be described in detail.

E_(seq)(C): Since different anatomical structures can be distinguishedbased on distinct temporal dynamics of the intensity-time curves, thisterm exploits the temporal information across the image sequence. Firsta metric is introduced that characterizes the distance between twodynamic signals to quantify the dissimilarity between different types oftemporal responses. Recall that I_(k)(x, y) denotes the intensity atpixel (x, y) in frame k. According to an embodiment of the invention,the intensity-time curve at each pixel (x, y) can be represented using acolumn vector:

I(x, y)=[I ₁(x, y), I ₂(x, y), . . . , I _(K)(x, y)]^(T).

Let θ denote the angle between two arbitrary vectors I₁ and I₂.According to an embodiment of the invention, a distance function betweenthe two vectors can be defined as

$\begin{matrix}\begin{matrix}{{{dis}\left( {I_{1},I_{2}} \right)} = {{\sin \frac{\theta}{2}}}} \\{{= \sqrt{\frac{1 - {\cos \; \theta}}{2}}},}\end{matrix} & (5)\end{matrix}$

where cos θ is given by the correlation coefficient between the meanremoved vectors Ĩ₁ and Ĩ₂ (the mean removed vectors are the vectorsafter subtracting the mean intensity for each particular vector), as thetemporal dynamics are of interest instead of the absolute intensity. Thecorrelation coefficient is given by

$\begin{matrix}{{{\cos \; \theta} = \frac{\langle{{\overset{\sim}{I}}_{1},{\overset{\sim}{I}}_{2}}\rangle}{{{\overset{\sim}{I}}_{1}}_{2}{{\overset{\sim}{I}}_{2}}_{2}}},} & (6)\end{matrix}$

with

denoting the inner product and ∥ ∥₂ the Euclidean L₂ norm.

According to an embodiment of the invention, it is assumed that thereare two regions in the image domain whose temporal signals closelyfollow two different dynamic profiles. Let Ω_(in)=inside(C) andΩ_(out)=outside(C) represent the regions inside and outside curve C,respectively. Let Ī_(prior) denote the prior known averageintensity-time curve for the object (inside curve C). In endocardium(epicardium) segmentation, Ī_(prior) can be estimated using the initialmask for the blood pool (myocardium).

To avoid meaningless noise incurred by regions exhibiting littletemporal dynamics, a spatially varying weight is assigned to each pixelthat lies between 0 and 1. The weight w(x, y) is defined as

$\begin{matrix}{{{w\left( {x,y} \right)} = \frac{Q\left( {x,y} \right)}{m\; {ax}_{({x,y})}{{Q\left( {x,y} \right)}}}},{{Q\left( {x,y} \right)} = {\sum\limits_{k = 1}^{K}\left\lbrack {{I_{k}\left( {x,y} \right)} - {I_{base}\left( {x,y} \right)}} \right\rbrack}},} & (7)\end{matrix}$

where I_(base)(x, y) is the baseline intensity for pixel (x, y),obtained by averaging the corresponding pixel intensities in the firstfew frames.

The energy term E_(seq)(C) is thus given by

$\begin{matrix}{{{E_{seq}(C)} = {{\int_{\Omega_{i\; n}}{{{w\left( {x,y} \right)}\begin{bmatrix}\begin{matrix}{{\alpha \; {{dis}^{2}\left( {{I\left( {x,y} \right)},{\overset{\_}{I}}_{i\; n}} \right)}} +} \\{{\beta \; {{dis}^{2}\left( {{I\left( {x,y} \right)},{\overset{\_}{I}}_{prior}} \right)}} +}\end{matrix} \\{\gamma \; {{dis}^{2}\left( {{\overset{\_}{I}}_{i\; n},{\overset{\_}{I}}_{prior}} \right)}}\end{bmatrix}}{x}{y}}} + {\int_{\Omega_{out}}{\left\lbrack {{w\left( {x,y} \right)}{{dis}^{2}\left( {{I\left( {x,y} \right)},{\overset{\_}{I}}_{out}} \right)}} \right\rbrack {x}{y}}}}},} & (8)\end{matrix}$

where Ī_(int) and Ī_(out) are the average intensity-time vectors forregion Ω_(in) and Ω_(out), respectively. α, β, and γ are positivescalars that lie in [0, 1], and satisfy α+β+γ=1. According to anembodiment of the invention, α may be chosen to be greater than 0.5. Thedistance functions are computed using EQS. (5) and (6). The integrals inEQ. (8) sum over the pixels in the respective regions, while thedistances in the integrands sum over the frames in the sequence.

E_(img)(C): This energy term exploits spatial (intra-image) informationin the reference frame I_(r). Assuming I_(r) is formed by two regions ofapproximately piecewise-constant intensities of distinct values, similarto EQ. (8), according to an embodiment of the invention, a region-basedenergy term E_(img)(C) can be defined as

$\begin{matrix}{{E_{img}(C)} = {{\int_{\Omega_{i\; n}}{\begin{bmatrix}\begin{matrix}{{\alpha \left( {{I_{r}\left( {x,y} \right)} - m_{i\; n}} \right)}^{2} +} \\{{\beta \left( {{I_{r}\left( {x,y} \right)} - m_{prior}} \right)}^{2} +}\end{matrix} \\{\gamma \left( {m_{i\; n} - m_{prior}} \right)}^{2}\end{bmatrix}{x}{y}}} + {\int_{\Omega_{out}}{\left( {{I_{r}\left( {x,y} \right)} - m_{out}} \right){x}{y}}}}} & (9)\end{matrix}$

where m_(in) and m_(out) are the average intensities inside Ω_(in) andΩ_(out), respectively; m_(prior) is the priori known average intensityfor the object (inside curve C). In endocardium (epicardium)segmentation, m_(prior) can be estimated using the initial mask for theblood pool (myocardium).

E_(int)(C): According to an embodiment of the invention, the contour maybe smoothed by minimizing its total arc length, using the internalenergy term:

E_(int)=∫_(C)ds   (10)

where ds represents the infinitesimal Euclidean arc length of thecontour.

E_(shape)(C): Because both endocardium and epicardium resemble ellipses,one can use an elliptical shape prior and penalize the squared distancebetween contour C and the ellipse estimated from C by direct leastsquare fitting of ellipses:

E _(shape) =∫ _(C) D ²(x, y)ds,   (11)

where D(x, y) is the distance between (x, y) and the estimated ellipse.

SEGMENTATION OF THE ENDOCARDIUM: An energy functional according to anembodiment of the invention for endocardium segmentation follows theformulation in EQ. (4). In particular, Ω_(in)=inside(C₁)=Ω_(a) andΩ_(out)=outside(C₁)=Ω_(b)∪Q_(c), as shown in FIG. 5.

Papillary muscles, which support the function of the valve, presentcertain issues. Papillary muscles are depicted in FIG. 3 by the circledarea 31 within an ROI 32. These muscles have the same intensity as themyocardium (heart wall) and are sometimes completely surrounded by(bright) blood, and sometimes blend into the wall. This can changeduring the heart cycle and a registration algorithm according to anembodiment of the invention should not be mislead, i.e. it shouldconsistently use the endocardial boundary. To address this issue,according to an embodiment of the invention, a weighted ellipse fittingmethod is used, in which each point on the contour was weightedaccording to its distance to the convex hull of the evolving contour C₁:the closer to the convex hull, the higher the weight.

SEGMENTATION OF THE EPICARDIUM: An energy functional according to anembodiment of the invention for epicardium segmentation has anadditional term that imposes anatomical constraints about the myocardiumthickness, that the distance between the endocardium and the epicardiumis relatively constant. In this case, C₁ is given from the results ofendocardium segmentation. Hence Ω_(in)=Ω_(b) and Ω_(out)=Ω_(c), as shownin FIG. 5.

Let d₁(x, y) denote the distance between point (x, y) and C₁. To takeadvantage of the prior knowledge that the distance between theendocardium and the epicardium is relatively constant, an energy termE_(con) is introduced as

E _(con)=∫_(C) ₂ f[d ₁(x, y)]ds,   (12)

where

$\begin{matrix}{{f\left\lbrack {d_{1}\left( {x,y} \right)} \right\rbrack} = \left\{ \begin{matrix}{0,{{{if}\mspace{14mu} d_{m\; i\; n}} \leq {d_{1}\left( {x,y} \right)} \leq d_{m\; {ax}}},} \\{{\infty,{{otherwise}.}}\mspace{149mu}}\end{matrix} \right.} & (13)\end{matrix}$

In EQ. (13), d_(min) and d_(max) are respectively the minimum and themaximum allowed myocardium thickness, which can be set according to theestimated myocardium thickness.

Results

FIG. 6 shows detected ROIs 61 overlaid on the reference frame for twocardiac MR perfusion sequences. In each sequence, the reference frameand the ROI were automatically determined by an algorithm according toan embodiment of the invention without any user input. FIG. 7 displaysthe segmentation results for the registered ROI sequence from FIG. 6,depicting the endocardial contour 71 and the epicardial contour 72.

System Implementations

It is to be understood that embodiments of the present invention can beimplemented in various forms of hardware, software, firmware, specialpurpose processes, or a combination thereof. In one embodiment, thepresent invention can be implemented in software as an applicationprogram tangible embodied on a computer readable program storage device.The application program can be uploaded to, and executed by, a machinecomprising any suitable architecture.

FIG. 8 is a block diagram of an exemplary computer system forimplementing a method for automatically generating magnetic resonance(MR) myocardial perfusion maps according to an embodiment of theinvention. Referring now to FIG. 8, a computer system 81 forimplementing the present invention can comprise, inter alia, a centralprocessing unit (CPU) 82, a memory 83 and an input/output (I/O)interface 84. The computer system 81 is generally coupled through theI/O interface 84 to a display 85 and various input devices 86 such as amouse and a keyboard. The support circuits can include circuits such ascache, power supplies, clock circuits, and a communication bus. Thememory 83 can include random access memory (RAM), read only memory(ROM), disk drive, tape drive, etc., or a combinations thereof. Thepresent invention can be implemented as a routine 87 that is stored inmemory 83 and executed by the CPU 82 to process the signal from thesignal source 88. As such, the computer system 81 is a general purposecomputer system that becomes a specific purpose computer system whenexecuting the routine 87 of the present invention.

The computer system 81 also includes an operating system and microinstruction code. The various processes and functions described hereincan either be part of the micro instruction code or part of theapplication program (or combination thereof) which is executed via theoperating system. In addition, various other peripheral devices can beconnected to the computer platform such as an additional data storagedevice and a printing device.

It is to be further understood that, because some of the constituentsystem components and method steps depicted in the accompanying figurescan be implemented in software, the actual connections between thesystems components (or the process steps) may differ depending upon themanner in which the present invention is programmed. Given the teachingsof the present invention provided herein, one of ordinary skill in therelated art will be able to contemplate these and similarimplementations or configurations of the present invention.

While the present invention has been described in detail with referenceto a preferred embodiment, those skilled in the art will appreciate thatvarious modifications and substitutions can be made thereto withoutdeparting from the spirit and scope of the invention as set forth in theappended claims.

1. A method of automatically generating a myocardial perfusion map froma sequence of magnetic resonance (MR) images, said method comprising thesteps of: determining a region of interest (ROI) in a reference frameselected from a time series of myocardial perfusion MR image slices;registering each image slice in said time series of slices to thereference frame to obtain a series of registered ROIs; and using theseries of registered ROIs to segment endo- and epi-cardial boundaries ofa myocardium in the ROI.
 2. The method of claim 1, wherein determining aregion of interest (ROI) in a reference frame comprises: cropping slicesin said series of slices according to a temporal intensity variation ateach pixel to obtain a global mask containing pixels exhibitingsignificant temporal intensity variation; selecting a range of candidateslices with superior contrast; detecting the LV blood pool and formingan LV mask in all candidate slices; selecting one candidate slice with adetected LV blood pool as the reference frame; detecting the rightventricle (RV) in the reference frame and forming an RV mask; estimatinga thickness of the myocardium based on the detected LV and RV masks; anddefining the ROI as a bounding box containing the LV, RV, andmyocardium.
 3. The method of claim 2, wherein detecting the LV bloodpool comprises, for each candidate slice, selecting an intensitythreshold for the blood, using said threshold to obtain a binary imageand removing small connected components from said binary image.
 4. Themethod of claim 3, wherein selecting one candidate slice as thereference frame comprises computing shape features for each remainingconnected component to select a least eccentric, most circular, and mostconvex connected component, wherein the selected connected componentamong all among all candidate slices corresponds to the LV blood pool.5. The method of claim 3, wherein said blood intensity threshold isselected through Gaussian clustering.
 6. The method of claim 4, whereinsaid shape features include an eccentricity computed as a ratio betweenminimum and maximum radii on the connected components, a convexitycomputed as a ratio between an area of the connected component and anarea of its convex hull, and a circularity of the convex hull defined as4π times the area of the convex hull over a square of a perimeter ofsaid convex hull.
 7. The method of claim 4, wherein when two connectedcomponents have similar shape features, a component with a larger areais selected.
 8. The method of claim 2, wherein detecting the RVcomprises finding a connected component in the binary image that isclose to the LV and has a large area.
 9. The method of claim 2, whereinsegmenting endo- and epi-cardial boundaries of a myocardium in the ROIcomprises: initializing the endo-cardium and epi-cardium boundaries fromthe LV mask and myocardial thickness estimate; detecting theendo-cardium by minimizing an energy functional; masking out aleft-ventricle (LV) blood pool inside of the endo-cardium from the imageseries, and detecting the epi-cardium by minimizing another energyfunctional, wherein each said energy functional incorporates correlationinformation among pixels in a same slice as well as temporalcorrelations across the slices in the time series, prior shapeinformation and anatomical constraints.
 10. The method of claim 9,wherein the energy functional for the endo-cardium isE(C₁)=E_(seq)(C₁)+λE_(img)(C₁)+μE_(int)(C₁)+κE_(shape)(C₁), wherein C₁is the endo-cardial boundary, λ, μ, and κ are parameters that controlweights of different terms, wherein${{E_{seq}(C)} = {{\int_{\Omega_{i\; n}}{{{w\left( {x,y} \right)}\begin{bmatrix}\begin{matrix}{{\alpha \; {{dis}^{2}\left( {{I\left( {x,y} \right)},{\overset{\_}{I}}_{i\; n}} \right)}} +} \\{{\beta \; {{dis}^{2}\left( {{I\left( {x,y} \right)},{\overset{\_}{I}}_{prior}} \right)}} +}\end{matrix} \\{\gamma \; {{dis}^{2}\left( {{\overset{\_}{I}}_{i\; n},{\overset{\_}{I}}_{prior}} \right)}}\end{bmatrix}}{x}{y}}} + {\int_{\Omega_{out}}{\left\lbrack {{w\left( {x,y} \right)}{{dis}^{2}\left( {{I\left( {x,y} \right)},{\overset{\_}{I}}_{out}} \right)}} \right\rbrack {x}{y}}}}},$wherein${{w\left( {x,y} \right)} = \frac{\sum\limits_{k = 1}^{K}\left\lbrack {{I_{k}\left( {x,y} \right)} - {I_{base}\left( {x,y} \right)}} \right\rbrack}{m\; {ax}_{({x,y})}{{\sum\limits_{k = 1}^{K}\left\lbrack {{I_{k}\left( {x,y} \right)} - {I_{base}\left( {x,y} \right)}} \right\rbrack}}}},$I_(k)(x, y) is the intensity of a pixel (x, y) in a kth slice of thetime series, I_(base)(x, y) is a baseline intensity for the pixel (x,y), I(x, y)=[I₁(x, y), I₂(x, y), . . . , I_(K)(x, y)]^(T) is aintensity-time curve for each pixel (x, y), Ī_(in) and Ī_(out) areaverage intensity-time vectors for region Ω_(in) and Ω_(out),respectively, Ω_(in) and Ω_(out) are the regions inside and outsideboundary C₁, respectively, Ī_(prior) is a prior known averageintensity-time curve inside boundary C₁, $\begin{matrix}{{{dis}^{2}\left( {I_{1},I_{2}} \right)} = {\frac{1}{2} - \frac{\langle{{\overset{\sim}{I}}_{1},{\overset{\sim}{I}}_{2}}\rangle}{2{{\overset{\sim}{I}}_{1}}_{2}{{\overset{\sim}{I}}_{2}}_{2}}}} \\{= \frac{{{{\overset{\sim}{I}}_{1}}_{2}{{\overset{\sim}{I}}_{2}}_{2}} - {\langle{{\overset{\sim}{I}}_{1},{\overset{\sim}{I}}_{2}}\rangle}}{2{{\overset{\sim}{I}}_{1}}_{2}{{\overset{\sim}{I}}_{2}}_{2}}}\end{matrix}$ where Ĩ is a mean removed vector, and α, β, and γ arepositive scalars that lie in [0, 1] and satisfy α+β+γ=1,${{E_{img}(C)} = {{\int_{\Omega_{i\; n}}{\begin{bmatrix}\begin{matrix}{{\alpha \left( {{I_{r}\left( {x,y} \right)} - m_{i\; n}} \right)}^{2} +} \\{{\beta \left( {{I_{r}\left( {x,y} \right)} - m_{prior}} \right)}^{2} +}\end{matrix} \\{\gamma \left( {m_{i\; n} - m_{prior}} \right)}^{2}\end{bmatrix}{x}{y}}} + {\int_{\Omega_{i\; n}}{\begin{pmatrix}{{I_{r}\left( {x,y} \right)} -} \\m_{out}\end{pmatrix}^{2}{x}{y}}}}},$ wherein I_(r) is the reference frame,m_(in) and m_(out) are average intensities inside Ω_(in) and Ω_(out),respectively and m_(prior) is a priori known average intensity insidecurve C₁, E_(int)=∫_(C) ₁ ds, an arc length of boundary C₁, andE_(shape)=∫_(C) ₁ D²(x, y)ds, wherein D(x, y) is a distance between apixel (x, y) in an elliptical prior shape and the estimated boundary C₁.11. The method of claim 10, further comprising weighting each pixel onthe endo-cardial boundary according to its distance to a convex hull ofendo-cardial boundary, wherein pixels closer to the convex hull have agreater weight.
 12. The method of claim 9, wherein the energy functionalfor the epi-cardium isE(C₂)=E_(seq)(C₂)+λE_(img)(C₂)+μE_(int)(C₂)+κE_(shape)(C₂)+E_(con)(C₂),wherein C₂ is the epi-cardial boundary, λ, μ, and κ are parameters thatcontrol weights of different terms, wherein${{E_{seq}(C)} = {{\int_{\Omega_{i\; n}}{{{w\left( {x,y} \right)}\begin{bmatrix}\begin{matrix}{{\alpha \; {{dis}^{2}\left( {{I\left( {x,y} \right)},{\overset{\_}{I}}_{i\; n}} \right)}} +} \\{{\beta \; {{dis}^{2}\left( {{I\left( {x,y} \right)},{\overset{\_}{I}}_{prior}} \right)}} +}\end{matrix} \\{\gamma \; {{dis}^{2}\left( {{\overset{\_}{I}}_{i\; n},{\overset{\_}{I}}_{prior}} \right)}}\end{bmatrix}}{x}{y}}} + {\int_{\Omega_{out}}{\left\lbrack {{w\left( {x,y} \right)}{{dis}^{2}\left( {{I\left( {x,y} \right)},{\overset{\_}{I}}_{out}} \right)}} \right\rbrack {x}{y}}}}},$wherein${{w\left( {x,y} \right)} = \frac{\sum\limits_{k = 1}^{K}\left\lbrack {{I_{k}\left( {x,y} \right)} - {I_{base}\left( {x,y} \right)}} \right\rbrack}{m\; {ax}_{({x,y})}{{\sum\limits_{k = 1}^{K}\left\lbrack {{I_{k}\left( {x,y} \right)} - {I_{base}\left( {x,y} \right)}} \right\rbrack}}}},$I_(k)(x, y) is the intensity of a pixel (x, y) in a kth slice of thetime series, I_(base)(x, y) is a baseline intensity for the pixel (x,y), I(x, y)=[I₁(x, y), I₂(x, y), . . . , I_(K)(x, y)]^(T) is aintensity-time curve for each pixel (x, y), Ī_(in) and Ī_(out) areaverage intensity-time vectors for region Ω_(in) and Ω_(out),respectively, Ω_(in) and Ω_(out) are regions inside and outside boundaryC₂, respectively, and Ω_(in) excludes a region inside the endo-cardiumboundary C₁, Ī_(prior) is a prior known average intensity-time curveinside boundary C₂, $\begin{matrix}{{{dis}^{2}\left( {I_{1},I_{2}} \right)} = {\frac{1}{2} - \frac{\langle{{\overset{\sim}{I}}_{1},{\overset{\sim}{I}}_{2}}\rangle}{2{{\overset{\sim}{I}}_{1}}_{2}{{\overset{\sim}{I}}_{2}}_{2}}}} \\{= \frac{{{{\overset{\sim}{I}}_{1}}_{2}{{\overset{\sim}{I}}_{2}}_{2}} - {\langle{{\overset{\sim}{I}}_{1},{\overset{\sim}{I}}_{2}}\rangle}}{2{{\overset{\sim}{I}}_{1}}_{2}{{\overset{\sim}{I}}_{2}}_{2}}}\end{matrix}$ where Ĩ is a mean removed vector, and α, β, and γ arepositive scalars that lie in [0, 1] and satisfy α+β+γ=1,${{E_{img}(C)} = {{\int_{\Omega_{i\; n}}{\begin{bmatrix}\begin{matrix}{{\alpha \left( {{I_{r}\left( {x,y} \right)} - m_{i\; n}} \right)}^{2} +} \\{{\beta \left( {{I_{r}\left( {x,y} \right)} - m_{prior}} \right)}^{2} +}\end{matrix} \\{\gamma \left( {m_{i\; n} - m_{prior}} \right)}^{2}\end{bmatrix}{x}{y}}} + {\int_{\Omega_{out}}{\begin{pmatrix}{{I_{r}\left( {x,y} \right)} -} \\m_{out}\end{pmatrix}^{2}{x}{y}}}}},$ wherein I_(r) is the reference frame,m_(in) and m_(out) are average intensities inside Ω_(in) and Ω_(out),respectively and m_(prior) is a priori known average intensity insidecurve C₁, E_(int)=∫_(C) ₂ ds, an arc length of boundary C₂,E_(shape)=∫_(C) ₂ D²(x, y)ds, wherein D(x, y) is a distance betweenpixel (x, y) and the estimated boundary C₂, and E_(con)=∫_(C) ₂ =f[d₁(x,y)]ds, wherein d₁(x, y) is a distance between a point (x, y) on C₂ andC₁, and${f\left\lbrack {d_{1}\left( {x,y} \right)} \right\rbrack} = \left\{ {\begin{matrix}{0,{{{if}\mspace{14mu} d_{m\; i\; n}} \leq {d_{1}\left( {x,y} \right)} \leq d_{m\; a\; x}}} \\{{\infty,{otherwise}}\mspace{155mu}}\end{matrix},} \right.$ wherein d_(min) and d_(max) are minimum andmaximum allowed myocardium thicknesses, respectively.
 13. A method ofautomatically generating a myocardial perfusion map from a sequence ofmagnetic resonance (MR) volumetric images, said method comprising thesteps of: providing a time series of registered myocardial perfusion MRimage slices having reference image with a masked left ventricle (LV);detecting the right ventricle (RV) in the reference frame and estimatinga thickness of the myocardium based on the detected LV and RV;initializing the endo-cardium and epi-cardium boundaries from the LVmask and myocardial thickness estimate; segmenting the endo-cardium byminimizing an energy functional; masking out a left-ventricle (LV) bloodpool inside of the endo-cardium from the image series, and segmentingthe epi-cardium by minimizing another energy functional, wherein eachsaid energy functional incorporates correlation information among pixelsin a same slice as well as temporal correlations across the slices inthe time series, prior shape information and anatomical constraints. 14.The method of claim 13, wherein providing a time series of registeredmyocardial perfusion MR image slices comprises providing a time seriesof myocardial perfusion MR image slices, determining a region ofinterest (ROI) in a reference frame selected from said time series, andregistering each image slice in said time series of slices to thereference frame to obtain a series of registered ROIs, wherein theseries of registered ROIs are used to segment endo- and epi-cardialboundaries of a myocardium in the ROI.
 15. A program storage devicereadable by a computer, tangibly embodying a program of instructionsexecutable by the computer to perform the method steps for automaticallygenerating a myocardial perfusion map from a sequence of magneticresonance (MR) images, said method comprising the steps of: determininga region of interest (ROI) in a reference frame selected from a timeseries of myocardial perfusion MR image slices; registering each imageslice in said time series of slices to the reference frame to obtain aseries of registered ROIs; and using the series of registered ROIs tosegment endo- and epi-cardial boundaries of a myocardium in the ROI. 16.The computer readable program storage device of claim 15, whereindetermining a region of interest (ROI) in a reference frame comprises:cropping slices in said series of slices according to a temporalintensity variation at each pixel to obtain a global mask containingpixels exhibiting significant temporal intensity variation; selecting arange of candidate slices with superior contrast; detecting the LV bloodpool and forming an LV mask in all candidate slices; selecting onecandidate slice with a detected LV blood pool as the reference frame;detecting the right ventricle (RV) in the reference frame and forming anRV mask; estimating a thickness of the myocardium based on the detectedLV and RV masks; and defining the ROI as a bounding box containing theLV, RV, and myocardium.
 17. The computer readable program storage deviceof claim 16, wherein detecting the LV blood pool comprises, for eachcandidate slice, selecting an intensity threshold for the blood, usingsaid threshold to obtain a binary image and removing small connectedcomponents from said binary image.
 18. The computer readable programstorage device of claim 17, wherein selecting one candidate slice as thereference frame comprises computing shape features for each remainingconnected component to select a least eccentric, most circular, and mostconvex connected component, wherein the selected connected componentamong all among all candidate slices corresponds to the LV blood pool.19. The computer readable program storage device of claim 17, whereinsaid blood intensity threshold is selected through Gaussian clustering.20. The computer readable program storage device of claim 18, whereinsaid shape features include an eccentricity computed as a ratio betweenminimum and maximum radii on the connected components, a convexitycomputed as a ratio between an area of the connected component and anarea of its convex hull, and a circularity of the convex hull defined as4π times the area of the convex hull over a square of a perimeter ofsaid convex hull.
 21. The computer readable program storage device ofclaim 18, wherein when two connected components have similar shapefeatures, a component with a larger area is selected.
 22. The computerreadable program storage device of claim 16, wherein detecting the RVcomprises finding a connected component in the binary image that isclose to the LV and has a large area.
 23. The computer readable programstorage device of claim 16, wherein segmenting endo- and epi-cardialboundaries of a myocardium in the ROI comprises: initializing theendo-cardium and epi-cardium boundaries from the LV mask and myocardialthickness estimate; detecting the endo-cardium by minimizing an energyfunctional; masking out a left-ventricle (LV) blood pool inside of theendo-cardium from the image series, and detecting the epi-cardium byminimizing another energy functional, wherein each said energyfunctional incorporates correlation information among pixels in a sameslice as well as temporal correlations across the slices in the timeseries, prior shape information and anatomical constraints.
 24. Thecomputer readable program storage device of claim 23, wherein the energyfunctional for the endo-cardium isE(C ₁)=E _(seq)(C ₁)+λE _(int)(C ₁)+μE_(int)(C ₁)+κE _(shape)(C ₁),wherein C₁ is the endo-cardial boundary, λ, μ, and κ and K areparameters that control weights of different terms, wherein${{E_{seq}(C)} = {{\int_{\Omega_{i\; n}}{{{w\left( {x,y} \right)}\begin{bmatrix}\begin{matrix}{{\alpha \; {{dis}^{2}\left( {{I\left( {x,y} \right)},{\overset{\_}{I}}_{\; {i\; n}}} \right)}} +} \\{{\beta \; {{dis}^{2}\left( {{I\left( {x,y} \right)},{\overset{\_}{I}}_{prior}} \right)}} +}\end{matrix} \\{\gamma \; {{dis}^{2}\left( {{\overset{\_}{I}}_{i\; n},{\overset{\_}{I}}_{prior}} \right)}}\end{bmatrix}}{x}{y}}} + {\int_{\Omega_{out}}{\left\lbrack {{w\left( {x,y} \right)}{{dis}^{2}\left( {{I\left( {x,y} \right)},{\overset{\_}{I}}_{out}} \right)}} \right\rbrack {x}{y}}}}},$wherein${{w\left( {x,y} \right)} = \frac{\sum\limits_{k = 1}^{K}\left\lbrack {{I_{k}\left( {x,y} \right)} - {I_{base}\left( {x,y} \right)}} \right\rbrack}{m\; {ax}_{({x,y})}{{\sum\limits_{k = 1}^{k}\left\lbrack {{I_{k}\left( {x,y} \right)} - {I_{base}\left( {x,y} \right)}} \right\rbrack}}}},$I_(k)(x, y) is the intensity of a pixel (x, y) in a kth slice of thetime series, I_(base)(x, y) is a baseline intensity for the pixel (x,y), I(x, y)=[I₁(x, y), I₂(x, y), . . . , I_(K)(x, y)]^(T) is aintensity-time curve for each pixel (x, y), Ī_(in) and Ī_(out) areaverage intensity-time vectors for region Ω_(in) and Ω_(out),respectively, Ω_(in) and Ω_(out) are the regions inside and outsideboundary C₁, respectively, Ī_(prior) is a prior known averageintensity-time curve inside boundary C₁, $\begin{matrix}{{{dis}^{2}\left( {I_{1},I_{2}} \right)} = {\frac{1}{2} - \frac{\langle{{\overset{\sim}{I}}_{1},{\overset{\sim}{I}}_{2}}\rangle}{2{{\overset{\sim}{I}}_{1}}_{2}{{\overset{\sim}{I}}_{2}}_{2}}}} \\{= \frac{{{{\overset{\sim}{I}}_{1}}_{2}{{\overset{\sim}{I}}_{2}}_{2}} - {\langle{{\overset{\sim}{I}}_{1},{\overset{\sim}{I}}_{2}}\rangle}}{2{{\overset{\sim}{I}}_{1}}_{2}{{\overset{\sim}{I}}_{2}}_{2}}}\end{matrix}$ where Ĩ is a mean removed vector, and α, β, and γ arepositive scalars that lie in [0, 1] and satisfy α+β+γ=1,${{E_{img}(C)} = {{\int_{\Omega_{i\; n}}{\begin{bmatrix}\begin{matrix}{{\alpha \left( {{I_{r}\left( {x,y} \right)} - m_{i\; n}} \right)}^{2} +} \\{{\beta \left( {{I_{r}\left( {x,y} \right)} - m_{prior}} \right)}^{2} +}\end{matrix} \\{\gamma \left( {m_{i\; n} - m_{prior}} \right)}^{2}\end{bmatrix}{x}{y}}} + {\int_{\Omega_{out}}{\begin{pmatrix}{{I_{r}\left( {x,y} \right)} -} \\m_{out}\end{pmatrix}^{2}{x}{y}}}}},$ wherein I_(r) is the reference frame,m_(in) and m_(out) are average intensities inside Ω_(in) and Ω_(out),respectively and m_(prior) is a priori known average intensity insidecurve C₁, E_(int)=∫_(C) ₁ ds, an arc length of boundary C₁, andE_(shape)=∫_(C) ₁ D²(x, y)ds, wherein D(x, y) is a distance between apixel (x, y) in an elliptical prior shape and the estimated boundary C₁.25. The computer readable program storage device of claim 24, the methodfurther comprising weighting each pixel on the endo-cardial boundaryaccording to its distance to a convex hull of endo-cardial boundary,wherein pixels closer to the convex hull have a greater weight.
 26. Thecomputer readable program storage device of claim 23, wherein the energyfunctional for the epi-cardium isE(C ₂)=E _(seq)(C ₂)+λE _(img)(C ₂)+μE _(int)(C ₂)+κE _(shape)(C ₂)+E_(con)(C ₂), wherein C₂ is the epi-cardial boundary, λ, μ, and κ areparameters that control weights of different terms, wherein${{E_{seq}(C)} = {{\int_{\Omega_{i\; n}}{{{w\left( {x,y} \right)}\begin{bmatrix}\begin{matrix}{{\alpha \; {{dis}^{2}\left( {{I\left( {x,y} \right)},{\overset{\_}{I}}_{i\; n}} \right)}} +} \\{{\beta \; {{dis}^{2}\left( {{I\left( {x,y} \right)},{\overset{\_}{I}}_{prior}} \right)}} +}\end{matrix} \\{\gamma \; {{dis}^{2}\left( {{\overset{\_}{I}}_{i\; n},{\overset{\_}{I}}_{prior}} \right)}}\end{bmatrix}}{x}{y}}} + {\int_{\Omega_{out}}{\left\lbrack {{w\left( {x,y} \right)}{{dis}^{2}\left( {{I\left( {x,y} \right)},{\overset{\_}{I}}_{out}} \right)}} \right\rbrack {x}{y}}}}},$wherein${{w\left( {x,y} \right)} = \frac{\sum\limits_{k = 1}^{K}\left\lbrack {{I_{k}\left( {x,y} \right)} - {I_{base}\left( {x,y} \right)}} \right\rbrack}{m\; {ax}_{({x,y})}{{\sum\limits_{k = 1}^{K}\left\lbrack {{I_{k}\left( {x,y} \right)} - {I_{base}\left( {x,y} \right)}} \right\rbrack}}}},$I_(k)(x, y) is the intensity of a pixel (x, y) in a kth slice of thetime series, I_(base)(x, y) is a baseline intensity for the pixel (x,y), I(x, y)=[I₁(x, y), I₂(x, y), . . . , I_(K)(x, y)]^(T) is aintensity-time curve for each pixel (x, y), Ī_(in) and Ī_(out) areaverage intensity-time vectors for region Ω_(in) and Ω_(out),respectively, Ω_(in) and Ω_(out) are regions inside and outside boundaryC₂, respectively, and Ω_(in) excludes a region inside the endo-cardiumboundary C₁, Ī_(prior) is a prior known average intensity-time curveinside boundary C₂, $\begin{matrix}{{{dis}^{2}\left( {I_{1},I_{2}} \right)} = {\frac{1}{2} - \frac{\langle{{\overset{\sim}{I}}_{1},{\overset{\sim}{I}}_{2}}\rangle}{2{{\overset{\sim}{I}}_{1}}_{2}{{\overset{\sim}{I}}_{2}}_{2}}}} \\{= \frac{{{{\overset{\sim}{I}}_{1}}_{2}{{\overset{\sim}{I}}_{2}}_{2}} - {\langle{{\overset{\sim}{I}}_{1},{\overset{\sim}{I}}_{2}}\rangle}}{2{{\overset{\sim}{I}}_{1}}_{2}{{\overset{\sim}{I}}_{2}}_{2}}}\end{matrix}$ where Ĩ is a mean removed vector, and α, β, and γ arepositive scalars that lie in [0, 1] and satisfy α+β+γ=1,${{E_{img}(C)} = {{\int_{\Omega_{i\; n}}{\begin{bmatrix}\begin{matrix}{{\alpha \left( {{I_{r}\left( {x,y} \right)} - m_{i\; n}} \right)}^{2} +} \\{{\beta \left( {{I_{r}\left( {x,y} \right)} - m_{prior}} \right)}^{2} +}\end{matrix} \\{\gamma \left( {m_{i\; n} - m_{prior}} \right)}^{2}\end{bmatrix}{x}{y}}} + {\int_{\Omega_{out}}{\begin{pmatrix}{{I_{r}\left( {x,y} \right)} -} \\m_{out}\end{pmatrix}^{2}{x}{y}}}}},$ wherein I_(r) is the reference frame,m_(in) and m_(out) are average intensities inside Ω_(in) and Ω_(out),respectively and m_(prior) is a priori known average intensity insidecurve C₁, E_(int)=∫_(C) ₂ ds, an arc length of boundary C₂,E_(shape)=∫_(C) ₂ D²(x, y)ds, wherein D(x, y) is a distance betweenpixel (x, y) and the estimated boundary C₂, and E_(con)=∫_(C) ₂ f[d₁(x,y)]ds, wherein d₁(x, y) is a distance between a point (x, y) on C₂ andC_(1,) and${f\left\lbrack {d_{1}\left( {x,y} \right)} \right\rbrack} = \left\{ {\begin{matrix}{0,{{{if}\mspace{14mu} d_{m\; i\; n}} \leq {d_{1}\left( {x,y} \right)} \leq d_{m\; {ax}}}} \\{{\infty,{otherwise}}\mspace{149mu}}\end{matrix},} \right.$ wherein d_(min) and d_(max) are minimum andmaximum allowed myocardium thicknesses, respectively.